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Abstract 

The  paper  demonstrates  how  multi-period  portfolio  optimization  problems 
can  be  efficiently  solved  as  multi-stage  stochastic  linear  programs.  A  scheme 
based  on  a  blending  of  classical  Benders  decomposition  techniques  and  a  special 
technique,  called  importance  sampling,  is  used  to  solve  this  general  class  of 
multi-stage  stochastic  linear  programs.  We  discuss  the  case  where  stochastic 
parameters  are  dependent  within  a  period  as  well  as  between  periods.  Initial 
computational  results  are  presented. 
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1.  Introduction 


Methods  of  Operations  Research,  especially  Mathematical  Programming  methods,  are 
receiving  broader  acceptance  in  the  financial  industry.  The  increasing  complexities 
and  inherent  uncertainties  in  financial  markets  have  lead  to  the  need  of  mathematical 
models  supporting  the  decision  making  process.  This  paper  addresses  the  portfolio 
selection  problem.  Since  Markowitz  (1959)  [20],  several  models  have  been  developed 
that  allow  one  to  determine  portfolios  with  the  highest  expected  returns  for  a  given 
level  of  risk.  His  model  (and  certain  closely  related  ones)  require  the  solution  of  a 
quadratic  program.  Other  approaches  model  the  stochastic  nature  of  the  problem 
directly  as  a  stochastic  program.  For  example,  Mulvey  (1987)  [21]  and  Mulvey  and 
Vladimirou  (1989)  [22]  [23]  formulate  asset  allocation  problems  as  a  stochastic  network 
problem. 

The  use  of  stochastic  programming  techniques  has  been  hampered  until  recently 
by  the  sheer  size  of  practical  problems  when  they  are  restated  as  deterministic  linear 
problems.  To  solve  them  it  was  necessary  that  the  number  of  scenarios  representing 
uncertainties  be  kept  small.  Most  models  developed  so  far  have  been  single-stage  or 
single- period  models,  that  is  to  say  to  the  case  where  the  decision  making  process  and 
the  future  events  (foresight)  are  restricted  to  a  single  time  period.  Only  few  attempts 
have  been  made  to  solve  practical  multi-stage  decision  making  models  whose  future 
events  are  spread  over  several  periods. 

Multi-stage  planning  problems  can  often  be  formulated  as  linear  programs  with 
a  dynamic  matrix  structure  which,  in  the  deterministic  case,  appear  in  a  staircase 
pattern  of  blocks  with  non-zero  submatrices.  These  blocks  correspond  to  and  are 
different  for  different  time  periods.  In  the  stochastic  case,  the  blocks  of  coefficients 
and  right  hand  sides  in  different  time  periods  are  functions  of  several  parameters 
whose  values  vary  stochastically  with  dependent  and  independent  distributions  which 
we  assume  to  be  known.  The  resulting  problem  is  a  multi-stage  stochastic  linear 
program.  Even  for  problems  with  a  small  number  of  stochastic  parameters  per  stage 
the  size  of  multi-stage  problems  when  expressed  in  equivalent  deterministic  form  can 
get  so  large  as  to  appear  intractable.  The  simplest  case  and  most  studied  is  that  with 
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two  stages.  Stochastic  linear  programs  were  first  introduced  by  Dantzig  (1955)  [4] 
and  Beale  (1955)  [1].  Since  then  it  has  been  studied  by  many  authors,  some  recent 
references  are  Birge  (1985)  [3],  Ermoliev  (1983)  [10],  Frauendorfer  (1988)  [12],  Higle 
and  Sen  (1989)  [14],  Kail  (1979)  [19],  Pereira  et  al.  (1989)  [25],  Rockafellar  and  Wets 
(1989)  [28],  Ruszczynski(1986)  [29],  and  Wets  (1984)  [31].  See  Ermoliev  and  Wets 
(1988)  [11]  for  a  survey  of  different  ways  proposed  to  solve  the  stochastic  programs. 

A  new  approach  based  on  Benders  decomposition  and  importance  sampling  was 
introduced  by  Dantzig  and  Glynn  (1990)  [5]  and  developed  jointly  by  them  and  In- 
fanger  (1990)  [17].  Our  approach  turned  out  to  be  very  powerful.  We  demonstrated 
its  power  by  solving  several  practical  large-scale  stochastic  linear  programs  witli  nu¬ 
merous  stochastic  parameters.  Infanger  (1991)  [18]  and  Dantzig  and  Infanger  (1991) 
[7]  report  on  computational  results  of  large-scale  problems  with  up  to  52  stocliastic 
parameters,  where  the  deterministic  equivalent  problem  if  attempted  to  express  it 
explicitly  would  have  had  several  billions  of  constraints.  These  problems  were  two- 
stage  problems  or  belonged  to  a  restricted  class  of  multi-stage  problems  which  could 
be  reexpressed  in  the  two-stage  framework. 

2.  The  Multi-Period  Asset  Allocation  Problem 

In  this  paper  we  formulate  a  class  of  multi-period  financial  asset  allocation  problems 
(Mulvey  and  Vladimirou  (1989)  [22])  and  show  how  they  can  be  solved  by  adaptations 
of  multi-stage  stochastic  linear  programs  methodology  and  software. 

At  the  initial  time  period  1  a  certain  amount  of  wealth  is  available  to  a  decision 
maker  in  assets  i  =  1, . . .  ,n  and  in  cash  which  we  index  as  asset  n  -t-  1.  We  denote 
Xi,i  =  1, . . .  ,Ti-l- 1  to  be  the  dollar  value  of  the  initially  available  assets.  The  decision 
maker  has  to  decide  each  period  how  to  rearrange  his  portfolio  to  achieve  best  return 
on  his  initial  investment  over  time.  We  consider  the  problem  in  discrete  time  and 
define  time  steps  t  =  1, ...  ,7’,  e.g.  by  months,  with  T  being  the  end  of  the  planning 
horizon. 

At  each  time  period  t  the  investor  can  either  hold  on  to  asset  ?,  buy  more,  or  .sell 
off  part  (or  all)  of  asset  i.  We  denote  y-  the  amount  sold  of  asset  i  in  period  t  and  by 
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I-  the  amount  of  asset  i  in  period  t  held  on  to.  Selling  means  decreasing  the  value 
x\  of  asset  i  and  increasing  the  value  of  cash,  iJi+i-  Also  the  investor  has  the  choice 
of  using  his  resulting  cash  to  buy  certain  amounts  of  assets  i.  The  amount  bought  in 
period  t  is  denoted  with  z\. 

Buying  and  selling  causes  transaction  costs  which  we  assume  to  be  proportional 
to  the  amount  of  dollar  value  of  asset  traded.  We  denote  by  lOOi/i  the  transaction 
costs  (expressed  as  a  percentage)  associated  with  buying  one  unit  of  i  and  with  lOO/i, 
the  transaction  costs  (expressed  as  a  percentage)  associated  with  selling  off  1  unit  of 
aisset  i.  Buying  1  unit  of  asset  i  requires  1  +i/,  units  of  cash  and  selling  1  unit  of  asset 
i  results  in  1  —  /i,  units  of  cash. 

Through  buying  and  selling  the  investor  can  restructure  his  portfolio  in  each  time 
period  t.  Once  this  f-th  stage  decision  is  made,  the  holdings  z  =  1,. . .  ,n  +  1  can 
be  calculated.  The  shares  in  the  portfolio  is  then  kept  constant  till  the  next  time 
period.  The  value  of  x\  is  affected  by  the  returns  on  the  market.  For  example  a 
portfolio  x\  at  time  t  changes  its  value  to  R\x\  where  R\  denotes  the  return  factors 
from  period  t  to  period  <  +  1 . 

At  time  t,  when  the  decision  on  rearranging  the  portfolio  has  to  be  made,  returns 
Rj,  for  i  =  l,...,n  are  not  known  to  the  decision  maker  with  certainty.  Only  the 
return  on  cash,  R^+i  is  assumed  known.  However,  we  assume  we  know  the  probability 
distributions  of  /?■.  The  problem  is  of  the  “wait-and-see”  type.  While  the  decision  at  t 
has  to  be  made  on  the  bcisis  of  distributions  of  future  returns  /?[,  for  z  =  1, . . . ,  n,  t  = 
l,...,r,  the  values  of  prior  returns  R\,  i  =  l,...,n,  t  =  l,...,t  —  1  have  already 
been  observed.  We  denote  with  /?*  =  R\,  for  t  =  1,. . .  ,n  the  n-dimensional  random 
vector  with  outcomes  r‘(W{),  uzj  €  Ht,  with  p"*  the  corresponding  probability  and  fl; 
the  set  of  all  possible  outcomes  in  t.  The  random  returns  R\  of  period  t  -  .n  mutually 
dependent  and  dependent  on  the  random  parameters  of  the  previous  period. 

After  the  last  period  T  no  decision  is  made.  Only  the  value  of  the  portfolio  is 
determined  by  adding  all  values  of  assets  including  the  laist  p^Tiod  returns.  We  call 
this  value  .  The  goal  of  the  decision  maker,  however,  is  to  maximize  Eu(v^),  the 
expected  utility  of  the  value  of  the  portfolio  after  period  T.  The  utility  function 
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u{v^)  describes  the  way  the  investor  views  risk.  If  u{v^)  is  linear,  it  describes  risk 
neutrality,  if  u(u^)  is  concave,  it  models  risk  averseness.  Nonlinear  utility  functions 
require  non  linear  programming  techniques  for  the  solution  of  the  problem.  Our 
methodology  is  not  restricted  to  linear  problems.  However,  for  the  sake  of  ease  and 
computational  speed  we  approximate  the  nonlinear  function  by  a  piecewise  linear 
function  with  sufficiently  large  number  of  linear  segments. 

In  the  model  presented  here  we  do  not  consider  shortselling  of  assets,  although 
this  feature  could  be  incorporated  easily.  We  also  do  not  consider  borrowing  of  cash, 
which  also  could  be  incorporated  easily.  The  holdings  of  assets,  as  well  as  the  amounts 
of  assets  sold  or  bought  have  to  be  positive.  In  general  there  are  also  lower  (x)  and 
upper  (i)  bounds  on  holdings  as  well  as  on  amounts  of  assets  to  be  sold  {y,y)  or  to 
be  bought  {z,z)  which  are  given  by  the  investor  and/or  by  the  market.  E.g.  a  certain 
asset  may  only  be  available  up  to  a  certain  amount  or  an  investor  wants  to  have  a 
certain  asset  with  at  least  a  certain  amount  of  dollar  value  in  the  portfolio.  Therefore 
in  general  we  formulate  x\  <  a:-  <  s-,  j/j  <  yj  <  y‘,  <  2'  <  z\,  where  x-  >  0,  y'  >  0, 
S.i  >  0,  given  for  f  =  1 , . . . , n  -f-  I,  ^  =  I, . . . ,  T. 

We  can  now  state  the  model: 


t  =  i,...,r. 

1  i  ■ 

=  l,...,n 

+ 

1,  rfx®  given: 

-r,  X. 

+ 

+ 

yl  - 

_ 

0,  ?:  =  !.. 

. . .  n 

-rl-,Wrr.\ 

+ 

x' 

•^n  +  I 

— 

0, 

u:irjxj 

+ 

T 

= 

0, 

max 

Eu{v^) 

x[  <  x[  < 

<  y', 

< 

y\,  ^<zl<zl 

i  =  1, . . .  ,n,  /  =  1. 

. T 

We  describe  correlation  between  asset  returns  using  a  factor  model.  Using  factors 
is  common  in  the  financial  industry  (e.g.  Perold  (1981)  [27]),  hence  historical  data  of 
various  factors  are  commercially  available.  The  idea  of  the  factor  model  is  to  K'late 

the  vector  of  asset  returns  /?'  =  (/fj, . . . , /?„)'  to  factors  =  (V] . \\Y .  While 

the  number  of  assets,  n  is  large,  e.g.  a  model  should  be  able  to  handle  about  500  to 
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3000  assets,  the  number  of  factors  h  is  comparatively  small.  Factor  models  used  in 
the  financial  industry  typically  involve  no  more  than  20  different  time  series  called 
factors.  The  factor  matrix  F{n  x  h)  relates  R‘  to  V*  : 

R‘  = 

The  coefficients  of  the  factor  matrix  are  estimated  using  regression  analyses  on  his¬ 
torical  data.  By  linear  transformations  of  historical  factors  the  transformed  factors 
can  always  be  determined  in  such  a  way  that  the  factors  are  orthogonal.  These 
factors  can  then  be  interpreted  as  independent  random  parameters  assumed  nor¬ 
mally  distributed  or  log  normally  distributed.  Using  the  factor  model  stochastically 
dependent  returns  can  be  generated  in  the  computer  by  using  these  stochastically 
independent  factors.  We  denote  the  random  factor  by  also  denoted  as  vj, 

with  corresponding  probability  p(v-),  where  pivj)  =  prob  (V-*  =  vj). 

We  also  consider  inter-period  dependency.  For  example  we  may  wish  to  have  a 
higher  probability  of  having  a  high  rate  of  return  in  period  t  if  it  was  high  in  period 
t  —  1  than  if  it  was  low  in  period  t  —  1.  We  can  model  this  inter-period  dependency 
cis  a  Markovian  type  process  applied  directly  on  the  factors: 

The  value  of  factor  i  in  period  t  is  the  sum  of  the  value  of  factor  i  in  the  previous 
period  f  -  1  plus  some  independent  random  variation  of  the  factor  in  t,  denoted  by 
7/‘.  The  Markovian  type  model  can  be  estimated  based  on  historical  data.  Instead 
of  having  an  additive  effect  as  above  we  may  prefer  to  have  a  multiplicative  effect 
by  applying  the  Markovian  process  directly  to  the  logs  of  the  factors.  We  haven’t 
explored  this  alternative. 

3.  Multi-Stage  Stochastic  Linear  Programs 

As  one  can  now  see  easily,  the  multi-period  asset  model  proposed  fits  exactly  into  the 
framework  of  a  general  class  of  multi-stage  stochastic  linear  programs  with  recourse. 
The  factor  model  for  generating  dependent  returns  and  the  Markovian  process  for 
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inter-period  dependency  define  a  special  class  of  dependencies  between  stochastic 
parameters  which  we  will  exploit  to  solve  the  problem.  Before  doing  so  we  state  the 
general  problem  and  the  methodology  we  have  developed  to  solve  it. 

The  multi-stage  stochastic  linear  program  can  be  formulated  as  follows: 

min  z  = 

CtX,  -h  E{c2X'^  . .  -f  EicTx'tp- ‘^))---) 

sft 

AiXi 

-B^xx  -I-  A2x'^ 


=  6^^ 


tjt  e  =  2,...,r 

The  problem  is  the  stochastic  extension  of  a  deterministic  dynamic  linear  program. 
While  the  first  stage  parameters  Ci,  Ai,  bi  are  known  to  the  planner  with  certainly, 
the  parameters  of  stages  2, ...  ,7’  are  assumed  known  only  by  their  distribution.  We 
assume  uncertainty  in  the  coefficients  of  the  transition  matrices  t  =  1 , . . . ,  T  and 
the  right  hand  sides  6"',  t  =  I,...,/*  and  assume  the  coefficients  of  the  technology 
matrices  At,  t  =  2, . . .  ,T  and  the  objective  function  coefficients  Cj,  t  =  2, ...  .T  to 
be  known  with  certainty.  The  goal  of  the  planner  is  to  minimize  the  expected  value 
of  present  and  future  costs. 

The  underlying  “wait-and-see”  decision  making-process  is  as  follows:  The  deci¬ 
sion  maker  makes  a  first  stage  decision  xi  before  observing  any  outcome  of  random 
parameters.  Then  he  waits  until  an  outcome  of  the  second  stage  random  parameters 
gets  realized.  The  second  stage  decision  then  is  made  based  on  the  knowledge  of  the 
realization  a>2  but  without  observing  any  outcome  of  random  parameters  of  stages 
2,  ...,T,  and  so  forth.  As  the  state  (the  actual  outcome)  is  carried  forward  to  tin- 
following  period,  the  decision  tree  grows  exponentially  with  the  number  of  stages.  W  e 
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consider  discrete  distributions  of  random  parameters  with  finite  number  of  outcomes, 
e.g.  Ut  €  ,  fit  =  { 1 , . . . ,  A' t } ,  t  =  1,..  .,T.  With  A't  being  the  number  of  scenarios 

in  period  t,  the  total  number  of  scenarios  for  ail  T  stages  is  A't.  The  number 
A't  is  expected  to  be  large,  as  it  is  computed  by  the  crossing  of  the  sets  of  possible 
outcomes  of  the  different  random  parameters  within  a  period.  E.g.  the  dimension  of 
the  random  vector  in  period  t  is  hj  and  17^  contains  kf  elements;  then  Kt  =  Ojii 
For  example,  in  the  asset  allocation  problem,  consider  the  case  of  20  factors,  modeled 
as  random  parameters  with  5  outcomes  each:  the  number  of  scenarios  per  period  is 

w  10'“*.  If  there  are  3  periods,  then  the  total  number  of  scenarios  grows  to  10^®. 
The  dimensions  of  an  equivalent  linear  program  of  an  asset  allocation  problem  with 
a  universe  of  about  500  assets  is  approximately  5  •  10^  rows  and  1.5  •  10^*  columns. 
It  is  of  course  impossible  to  write  down  this  linear  program  explicitly. 

It  is  clear  that  the  multi-period  asset  allocation  problem  defined  above  is  a  special 
case  of  the  multi-stage  stochastic  linear  program.  The  correspondence  is  as  follows: 
the  vector  Xt  now  denotes  the  vector  of  all  decision  variables  (holdings,  amount  to 
be  bought  and  to  be  sold)  in  period  t.  Uncertainty  occurs  only  in  the  transition 
matrices  5,  which  contain  in  their  diagonal  the  return  factors  A-.  The  right  hand 
sides  62, . . . ,  are  zero,  as  well  as  the  objective  function  coefficients  C2, . . . ,  cj-i-  We 
now  describe  the  techniques  we  have  developed  to  solve  the  multi-stage  program. 

4.  Benders  Decomposition 

A  description  of  how  Benders  (1962)  [2]  Decomposition  Algorithm  can  be  applied  to 
solve  stochastic  linear  programs  can  be  found  in  Van  Slyke  and  Wets  (1969)  [30],  Birge 
(1985)  [3].  Using  Benders  decomposition  we  decompose  the  problem  into  subproblems 
of  different  stages  t.  In  the  most  general  case  where  there  is  a  dependency  of  stochastic 
parameters  between  stages  the  number  of  subproblems  is  equal  to  the  number  of 
scenarios  in  each  stage  t.  To  distinguish  one  subproblem  from  another,  each  is  indexed 
with  Ut, . . .  ,uj2,  where  u>t  is  the  random  event  in  stage  t  and  Wi-i, . . .  ,uj2  is  the  path 
of  previous  events  which  gave  rise  to  the  particular  subproblems  in  stage  t. 

For  expository  purposes,  we  assume  initially  the  random  events  that  happen  in 
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one  stage  are  independent  of  those  that  happen  in  the  next  stage.  For  example,  when 
the  probability  of  having  a  high  rate  of  return  in  period  t  is  the  same  for  all  values  of 
rate  of  return  in  period  t  —  l.  In  the  independent  case  scenarios  G  Ht+i  in  period 
t  +  1  are  identical  for  each  scenario  Ut  €  in  period  t.  The  history  is  only  carried 
forward  through  optimal  decisions  from  previous  periods.  In  the  special 

class  of  Markovian  dependency  which  we  described  earlier,  ■> 

where  tt  represents  a  matrix  of  random  parameters  independent  of  those  in  period 
t  -  1. 

The  idea  of  using  Benders  decomposition  is  to  express  in  each  stage  t,  t  = 
1, . . . ,  T—l  and  scenario  wt  the  expected  future  costs  (the  impact  of  stages  <  +  l, . . . ,  7’) 
by  a  scalar  6t  and  “cuts”,  necessary  conditions  for  feasibility  and  optimality  which  are 
expressed  only  in  terms  of  the  stage  t  decision  variables  Xt  and  Ot.  Cuts  are  initially 
absent  and  then  sequentially  added  to  the  stage  t  problems.  Each  scenario  subprob¬ 
lem  uji  in  stage  t  collects  the  information  about  expected  future  costs  by  means  of 
the  cuts. 

The  relation  between  the  stages  and  scenarios  in  the  decomposed  multi-stage 
problem  is  summarized  as  follows: 

Stage  1  problem: 


min  Zj  = 

Cili 

+ 

Oi 

sjt 

TTi  : 

AiXi 

= 

bx 

pi'  : 

H 

1 

+ 

> 

II 

I, 

> 

Ox 

> 

0 

=  2,...,T- 

1,  problem: 

min  z"'  = 

CtXt‘ 

+ 

sjt 

: 

AtxT 

= 

bT  + 

(j(  X( 

+ 

eT 

> 

II 

xf 

1 

9T 

> 

0 
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Stage  T  problem: 


•  Uf'P 

mm  Zj  = 

U/T 

CtXt 

s/t 

_"T  . 

Ti'j' 

Ajx'tp' 

=  651^  + 

IV 

0 

min  2i  represents  the  optimal  objective  function  value  in  the  first  stage,  ii,  6i 
represent  the  optimal  solution,  the  vector  tti  denotes  the  optimal  dual  prices  associ¬ 
ated  to  the  original  stage  1  constraints,  and  the  scalars  pj*  are  the  optimal  dual  prices 
associated  to  the  cuts,  which  have  been  added  so  far  in  iterations  l\  = 

The  optimal  objective  function  values  min  2"'  =  min  and  the  optimal 

dual  prices  tt"'  =  7r“'(xf_i)  associated  to  the  original  stage  t  constraints  in  stages 
f,  t  =  2, . . .  ,T  and  the  optimal  dual  prices  associated  to  the  cuts 

in  stages  <,  t  =  2, . . .  ,T  —  1  are  all  dependent  upon  the  optimal  solution  passed 
as  input  from  the  previous  stages  t  —  l.  According  to  the  scenario  development  in  the 
previous  stages  an  optimal  solution  is  actually  indexed  by  the  scenario  outcomes 
of  all  previous  stages  and  is  therefore  denoted  as  For  exposition, 

we  suppress  the  scenario  history  and  present  the  optimal  solution  of  subproblems  in 
stage  t,  scenario  as  a  function  of  the  input  it_i. 

We  compute  the  expected  future  costs  as  zt+i  =  the  right  hand  sides 

of  the  cuts  as  gl'  =  coefficients  of  the 

cuts  as  G'/  =  where  p^^  =  0,  G'^'^  =  0,  and  g'^'  =  0. 

A  subproblem  in  stage  t  and  in  scenario  uJt  interacts  with  its  predecessors  and 
descendants  by  passing  forward  optimal  solutions  and  backwards  cuts.  Benders  de¬ 
composition  splits  the  multi-stage  problem  into  a  series  of  two-stage  relations  which 
are  overall  connected  by  a  nesting  scheme.  We  call  the  stage  f,  scenario  u;^  problem 
the  current  master  problem.  It  receives  from  its  ancestor  in  period  <  —  1  a  solution 
it_i.  The  current  scenario  is  determined  by  the  outcome  u>t  of  the  random  param¬ 
eters  in  stage  t  which  are  reflected  in  the  right  hand  side  6"*  T  Bf'liXt-i.  As  stated 
above,  Z(_i  hais  a  history.  The  history  has  to  be  considered  when  nesting  several 
stages.  Given  and  subject  to  Xt-i  we  solve  the  stage  t  problem  in  scenario  Ut  and 
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pass  the  obtained  solution  x"‘  to  the  descendant  problems.  By  so'^  ing  all  problems 
€  fit+i  (referred  to  as  the  universe  case)  we  compute  the  expected  value  of  the 
descendant  stage  costs  zt+i  =  Eu,+\Z^^\'  and  the  coefficients  Gt  = 
and  the  right  hand  side  gt  =  9t+i)  of  a  cut.  The 

cut  is  added  to  the  current  master  problem  (stage  t,  scenario  problem)  and  by 
solving  the  problem  again  another  trial  solution  is  obtained. 

The  optimal  solution  the  of  current  master  problem  in  stage  t,  scenario  lOt  gives 
a  lower  bound,  and  the  expected  cost  of  the  trial  solution  gives  an  upper  bound  of 
the  expected  costs  of  all  sceiiarios  descendant  from  the  stage  t  scenario  Ut-  If  lower 
bound  and  upper  bound  are  sufficiently  close,  the  current  master  problem  is  said 
to  represent  the  future  expected  cost  and  contains  (by  means  of  a  sufficient  number 
of  cuts)  all  the  information  needed  from  future  scenarios.  In  this  case  w’e  say  the 
current  master  is  balanced  with  its  descendant  problems. 

Note  that  the  current  master  problem  represents  the  expected  future  costs  only 
subject  to  the  trial  solution  it_i  which  was  passed  from  its  ancestor  and  subject  to  the 
current  scenario  Wj.  Note  also  that  we  have  implicitly  assumed  that  the  descendant 
problems  in  stage  t  + 1  are  also  balanced  with  their  descendant  problems  in  stage  t  +  2 
by  means  of  having  collected  a  sufficient  number  of  cuts  to  represent  the  '^xpect 'd 
costs  of  descendant  scenarios  for  <+2  on,  and  so  forth.  However,  note  that  the  solution 
of  the  current  stage  t  scenario  Ut  problem  gives  a  lower  bound  of  the  expected  costs 
of  all  scenarios  descendant  from  the  stage  t  scenario  coj  problem  regardless  of  having 
collected  a  sufficient  number  of  cuts.  We  shall  exploit  this  fact. 

Two  properties  of  cuts  are  crucial  for  the  solution  procedure; 

1.  In  the  case  of  independence  of  stochastic  parameters  between  stages: 
The  cuts  derived  from  any  trial  solution  are  valid  cuts  for  all  subproblems  u,’,  G 
Ot.  E.g.  the  cut:  6t  >  Xt  +  i a  constraint  who.s<> 

coefficients  don’t  depend  on  x,,  hence  is  valid  for  all  values  of  X(.  To  see  this,  noli- 

optimal  dual  prices  that  do  depend  on  .r;  for  optimality  but 
they  remain  dual  feasible  independent  of  the  values  of  the  right  hand  side  as  a 
function  of  X(.  The  validity  of  the  cuts  depends  only  on  the  dual  feasibility  of  the 


11 


^r+V  •  represents  an  outer  linearization  of  the  future  expected  cost  function  2(+i(X(), 
evaluated  at  if  Different  scenarios  Wj  are  in  stage  t  are  distinguished  by  different 
right  hand  sides  of  the  original  stage  t  constraints,  e.g.  AtXt  =  The 

set  of  cuts  —G[‘xt  +  0t  >  gi'ih  =  represent  an  outer  linearization  of  the 

expected  future  costs  independent  of  scenarios  ujt  G  ^t-  The  outer  linearization 
defined  by  the  set  of  cuts  equals  the  expected  future  cost  function,  if  (xj)  =  6^^ 

where  Ot  is  the  value  of  dt  corresponding  to  the  solution  it  of  any  stage  t  problem, 
ff  Ez‘^^y(i'f')  =  ut  6  Di,  then  a  sufficient  number  of  necessary  cuts  have  been 
generated  to  represent  the  expected  future  costs  for  all  solutions  x"'  of  scenarios 
ujt  G  Dj  in  stage  t  and  we  say  stage  t  is  balanced  with  stage  <  +  1. 

2.  In  the  case  of  dependency  of  stochastic  parameters  between  stages: 
Cuts  now  depend  on  scenario  Ut  in  period  t.  Sharing  of  cuts  between  different  scenario 
subproblems  G  Dj  is  no  longer  directly  possible.  However,  for  additive  dependency, 
(e.g.  Markovian  type  dependency)  cuts  can  be  easily  adjusted  to  different  scenarios. 
(See  Pereira  and  Pinto  (1989)  [26]  for  additive  dependent  right  hand  sides.)  For 
example  in  the  case  of  the  Markovian  type  dependency  which  we  introduced  in  the 
multi-period  asset  allocation  problem  -f  Here  t,  represents  a 

matrix  whose  elements  are  functions  of  random  parameters  which  are  independent  of 
the  period  t  — 1  random  parameters.  (The  elements  of  e  are  the  independent  part  of  the 
random  returns  and  are  generated  by  the  product  F  rji  where  rjt  is  the  random  change 
in  Vt  that  generated  C(.)  In  the  case  of  the  additive  dependency  a  cut  in  stage  t  and 
scenario  wt  has  the  form:  6t  >  [(£'w,+, -b  E^.+i  e"“^']x£ -f  . 

It  can  be  easily  seen  that  the  coefficients  of  the  cut  consist  of  a  part  independent  of 
scenarios  Ut  and  a  dependent  part.  The  cut  can  be  adjusted  to  different  scenarios 
u>t  £  Dj  by  adding  the  scenario  dependent  part  {E^,^jw'^^\')B'^‘  according  to  scenario 
ujf  This  requires  storing  of  the  expected  value  of  the  dual  variables  • 

Taking  advantage  of  the  above  stated  properties  we  actually  only  need  to  store 
one  subproblem  per  stage  t.  For  different  scenarios  Ut  and  different  solutions  X(_i 
passed  from  the  previous  stage  we  determine  the  right  hand  side  accordingly.  The 
cuts  are  valid  for  all  scenarios  Ut  G  Dj  in  the  case  of  independence  of  the  stochastic 
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parameters  between  stages  or  are  adjusted  in  the  gradient  according  to  the  actual 
scenario  Ut  in  case  of  Markovian  type  dependency  between  stages.  Therefore  it  is 
easily  possible  to  generate  any  u>t  subproblem.  Future  information  is  represented  in 
the  cuts  which  have  been  generated  so  far  and  can  be  efficiently  used  in  any  scenario 
u>t  G  fit  independently  from  which  scenario  originated  it. 

5.  Multidimensional  Integration 

The  computation  of  the  expected  future  costs  Zt+i,  the  coefficients  Gt  and  the  right 
hand  side  gt  of  the  cuts  requires  the  computation  of  multiple  integrals  or  multiple 
sums.  The  expected  value  of  the  second  stage  costs  in  period  <  +  1  (we  suppress  the 
index  t  for  this  discussion),  e.g.  z  =  Ez'^  =  E(C)  is  an  expectation  of  functions 
C(u“'),u;  G  n,  where  C(u")  is  obtained  by  solving  a  linear  program.  V  (in  general)  is 
a  h-dimensional  random  vector  parameter,  e.g.  V  =  ( Vi, . . . ,  V/,),  with  outcomes  v'^  = 
(i>i,. . .  ,U/,)'^.  For  example  V  represents  the  value  of  the  i-th  factor  u"  the  observed 
random  outcome.  The  vector  u"  is  also  denoted  by  u,  and  p(u'*')  alias  p(u)  denotes  the 
corresponding  probability,  is  the  set  of  all  possible  random  events  and  is  constructed 
by  crossing  the  sets  of  outcomes  =  fii  x  Q2  x  •  •  •  x  With  P  being  the  probability 
measure  under  the  assumption  of  independence  the  integral  E  C(V')  =  /  C{v'^)P{duj) 
takes  the  form  of  a  multiple  integral  E  C{V)  =  f  ■  •  ■  f  C(v)p(v)dvi  ■  ■  ■  dr/,,  or.  in  case 
of  discrete  distributions,  the  form  of  a  multiple  sum  E  C(V)  =  Ev,  ’ '  ’  E,.,  C(r)p(r), 
where  p(r)  =  pi(ri)  •  •  •p/.(r/,). 

The  number  of  terms  in  the  multiple  sum  computation  gets  astronomically  large 
and  therefore  the  evaluations  of  multiple  sums  by  direct  summation  is  not  practical. 
This  is  especially  true  because  function  evaluations  are  computationally  expensive 
since  the  evaluation  of  each  term  in  the  multiple  sum  requires  the  solution  of  a  linear 
program.  In  the  following  we  discuss  a  scheme  for  estimating  the  exi)ected  vahu's 
with  a  sufficiently  low  estimation  error  without  having  to  evaluate  each  t('rm. 
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6.  Importance  Sampling 


Monte  Carlo  Methods  are  recommended  to  compute  multiple  integrals  or  multiple 
sums  for  higher  /i-dimensional  sample  spaces  (Davis  and  Rabinowitz  (1984)  [9]  ,  Glynn 
and  Iglehart  (1989)  [13]).  Suppose  =  (^(u")  are  independent  random  variates 
of  n",  ui  =  l,...,n  with  expectation  z,  where  n  is  the  sample  size.  An  unbiased 
estimator  of  z  with  variance  =  var{C{V))  is 


i  =  (i/n)y;c“. 


u;=l 

Note  that  the  standard  error  decreases  with  ®  and  the  convergence  rate  of  2  to  r  is 
independent  of  the  dimension  of  the  sample  space  h.  We  rewrite  2  = 
as 


E 


C{v‘^)p{v'^)q{v'^) 


by  introducing  a  new  probability  mass  function  ^(v")  and  we  obtain  a  new  estimator 

1  "  C{v'^)p{v‘^) 


=  rE 


n 


u/=l 


q{v-’) 


by  sampling  from  9(v").  The  variance  of  2  is  given  by 


Choosing  q’iv'^)  =  ....  would  lead  to  var{z)  =  0,  which  means  one  could 

get  a  perfect  estimate  of  the  multiple  sum  from  only  one  estimation.  Practically, 
however,  this  is  useless  since  to  compute  ^(u*^)  we  have  to  know  2  =  C"^p(*’'^)i 

which  is  what  we  are  trying  to  compute  in  the  first  place. 

The  result,  however,  helps  to  derive  a  heuristic  for  choosing  q.  It  should  be 
proportional  to  the  product  C'(v")p(u‘^)  and  should  have  a  form  that  can  be  integrated 
easily.  Thus  a  function  r(v‘^)  w  C{v‘^)  is  sought,  which  can  be  integrated  with  less 
effort  than  C(u").  Additive  and  multiplicative  (in  the  components  of  the  stochastic 
vector  v)  approximation  functions  and  combinations  of  these  are  potential  candidates 
for  our  approximations.  Especially  for  financial  investment  problems,  we  have  been 
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getting  good  results  using  C(V)  fs  C,(K)-  We  compute  q  as 


C{v'^)p{v'^) 

I2«=l  I3uf€n, 


In  this  case  one  has  to  compute  only  h  1-dimensional  sums  instead  of  1  /i-dimensional 
sum.  The  variance  reduction  depends  on  how  well  the  approximation  function  fits 
the  original  cost  function.  If  the  original  cost  function  has  the  property  of  additivity 
(separability)  the  multiple  sum  can  be  computed  exactly  by  h  1 -dimensional  sums.  If 
the  additive  model  is  a  bad  approximation  of  the  cost  function  the  only  “price”  that 
has  to  be  paid  is  increasing  the  sample  size.  If  the  observed  variance  is  too  high  using 
a  starting  sample  size,  the  sample  size  is  adjusted  higher.  Actually  we  use  a  variant  of 
the  additive  approximation  function.  By  introducing  C(t),  the  costs  of  a  base  case, 
we  make  the  model  more  sensitive  to  the  impact  of  the  stochastic  parameters  v. 

h 


r(r)  =  c{r)  +  Y.  r.lK),  r.(K)  =  c(t„  . . 


A-1 5  A'+i  1  •  •  • ' 


1=1 

We  denote  this  as  a  marginal  cost  model,  r  can  be  any  arbitrary  chosen  point  of  the 
set  of  values  u,,  ?  =  1, . . .  ,/i.  For  example  we  choose  r,  as  that  outcome  of  V’,  which 
leads  to  the  lowest  costs,  ceteris  paribus. 

Summarizing,  the  importance  sampling  scheme  has  two  phases:  the  preparation 
phase  and  the  sample  phase.  In  the  preparation  pheise  we  explore  the  cost  function 
C(V)  at  the  margins  to  compute  the  additive  approximation  function  r(V’).  For  this 
process  n^rep  =  1  +  —  1)  subproblems  have  to  be  solved.  Using  FfV')  we 

compute  the  approximate  importance  density 


_ r(u“')p(u“') _ 

c(r)  +  !:!•=,  Z!u;6n, 


Next  we  sample  n  scenarios  from  the  importance  density  and,  in  the  sample  phase, 
solve  n  linear  programs  to  compute  the  estimation  of  r  using  the  Monte  Carlo  esti¬ 
mator.  We  compute  the  gradient  G  and  the  right  hand  side  q  of  the  cut  using  the 
same  sample  points  at  hand  from  the  expected  cost  calculation.  See  Infangi'r  (19hl) 
[18]  for  the  computation  of  the  cuts  and  details  of  the  estimation  process. 
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7.  The  Algorithm 

By  solving  a  sample  of  subproblems  Wj+i  according  to  the  importance  sampling  scheme 
we  compute  estimates  of  the  expected  future  costs  the  gradients  G[‘  and  the 
right  hand  sides  g'/  of  the  cuts  in  each  stage  t  and  scenario  ujf  The  objective  function 
value  of  the  solution  of  each  stage  t,  scenario  Wt  subproblem  gives  a  valid  lower-bound 
estimate  of  the  expected  costs  2"'  =  Cti"'  -f-  subject  to  scenario  u;*  and  subject  to 
it-i,  the  (optimal)  solution  passed  forward  from  the  previous  stage.  The  obtained 
lower-bound  estimate  is  the  tightest  lower  bound  that  can  be  generated,  if  in  stage 
t-|- 1  a  sufficient  number  of  cuts  have  been  added  to  represent  the  expected  future  costs 
with  respect  to  stage  #  -t-  1  for  all  scenarios  u>(+i  G  Ot+i  and  is  a  weaker  lower-bound 
estimate  if  there  is  not  a  sufficient  number  of  cuts. 

We  are  especially  interested  in  the  lower-bound  estimate  of  the  first  stage  costs 
which  we  obtain  by  solving  the  first  stage  problem.  If  the  first  stage  problem  is 
balanced  with  the  second  stage,  that  is,  if  the  cuts  added  so  far  to  the  first  stage 
problem  fully  represent  the  expected  second  stage  costs,  and  if  the  second  stage  is 
balanced  with  the  third  stage  for  all  scenarios  u>2  €  ^2  and  all  values  of  ii,  passed 
to  it  from  the  first  stage,  and  so  forth  till  stage  T  —  1,  then  the  solution  of  the  first 
stage  problem  is  the  optimum  solution  of  the  multi-stage  stochastic  linear  program. 
In  this  case  the  lower  bound  estimate  of  Zi  takes  on  the  value  of  the  total  expected 
costs  of  the  multi-stage  problem. 

To  obtain  an  upper  bound  of  the  total  expected  costs  of  the  multi-stage  problem, 
we  evaluate  the  expected  costs  of  the  current  first  stage  trial  solution  f  j.  This  can  be 
accomplished  by  sampling  paths  from  stages  2, . . . ,  T.  For  a  reference,  see  Pereira  and 
Pinto  (1989)  [26].  To  efficiently  sample  a  small  number  of  paths  to  obtain  an  accurate 
estimate  of  the  expected  costs  associated  with  ij,  we  also  use  importance  sampling. 
We  define  a  path  s"  =  (xj ,  X2, . . . ,  x?-)",  u;  €  where  fl  =  {^2  x  fia  x  •  •  •  x  fir}, 
as  a  sequence  of  optimal  solutions  x'^'  of  stage  t  scenario  uit  problems,  t  =  2, . . . ,  T 
and  X]  being  the  first  stage  trial  solution.  A  path  is  computed  by  observing  the 
“wait-and-see”  requirements;  We  pass  Xi  to  the  second  stage  and  solve  the  second 
stage  problem  for  scenario  u>2  and  obtain  the  optimal  solution  x^.  Next  we  pass  the 


16 


obtained  second  stage  solution  to  the  third  stage  and  solve  the  third  stage  problem 
for  scenario  wa  to  obtain  x^.  We  continue  in  this  way  until  we  obtain  x'^  in  stage 
T.  Note  that  when  solving  the  stage  t  problem  no  future  outcomes  Wj+i, . . .  ,ujt  are 
used.  All  future  information  at  each  stage  is  solely  represented  by  means  of  the  cuts 
added  in  stage  t  so  far.  The  costs  of  a  path  i",  0(3")  is  given  by  C{s'^)  =  ■ 

The  expected  value  of  the  costs  of  all  paths  i",  a;  €  fi,  Es'^  gives  an  upper  bound  to 
the  costs  of  a  trial  solution  Xi. 

We  sample  paths  by  applying  the  importance  sampling  scheme  to  the  dimensional 
space  of  size  YlJ=2^t  random  parameters  V^|,  it  =  t  =  2,...,7'. 

For  sampling  paths  the  importance  density  ^(V')  is  computed  based  on  the  additive 
marginal  approximation  function  analogous  to  the  way  it  was  defined  earlier: 


r(K)  =  C{r)  +  5:i:C(r.., . r. 

,ii-l »  ^t,ii )  +  l '  •  •  •  '  Tr.hr)  - 


i,  =  l 


where  V  =  •  •  • ,  ,  V;^, . . . ,  and  r  =  (t/,  . . .  ,7)}, ,  r^, . . . ,  .  Sampling 

paths  u;  €  n  according  to  this  importance  sampling  scheme  we  obtain  an  equal 
number  of  sample  points  u>t  €  Clt  in  stages  t  =  2,...,7’.  At  these  sample  points  we 
define  the  current  stage  t  scenario  ut  subproblems  and  generate  cuts  to  be  added  at 
stages  <  =  1,...,T—  Iby  employing  importance  sampling  as  described  above  for 


The  overall  procedure  works  as  follows:  Solving  the  stage  1  problem  in  iteration 
1  we  obtain  a  trial  solution  xi  and  a  lower  bound  estimate  of  the  expected  costs  cj. 
Now  we  employ  the  path  sampling  procedure  to  obtain  an  upper  bound  estimate  of 
the  expected  costs  Zj.  If  the  upper  bound  estimate  and  the  lower  bound  estimate 
are  within  a  given  optimality  tolerance,  we  call  the  first  stage  solution  the  optimal 
solution  of  the  multi-stage  problem,  and  quit.  Otherwise  we  generate  cuts  in  stages 
1, . . . ,  r  —  1.  The  path  sampling  procedure  used  for  the  upper  bound  estimate  has 
produced  sample  points  u>t  €  in  stages  t  =  2,  ...,T  with  corresponding  ancestor 
solutions  Xi  and  x"'  in  stages  t  =  2, . . . ,  T  —  1  to  be  passed  to  the  current  stage 
t  scenario  Ut  problem.  Starting  at  stage  T  —  1  and  moving  backwards  till  stage  1 
we  take  each  sample  problem  Ut  in  stage  t  and  finally  the  stage  I  problem  as  the 
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current  master  problem  and  compute  cuts  by  sampling  again  Wj+i  €  Hj+i  descendant 
problems  until  each  scenario  problem  Ut  in  stage  t  is  balanced  with  stage  <  +  1  with 
regard  to  ancestor  solutions  which  have  been  passed  from  stage  <  —  1.  Arriving 
at  stage  1  we  obtain  a  new  solution  xi  and  a  new  lower  bound  estimate.  We  continue 
as  defined  above  by  sampling  new  paths  for  the  upper  bound  estimate.  Finally,  after 
a  finite  number  of  iterations,  upper  and  lower  bound  estimates  will  be  sufficiently 
close.  Upper  and  lower  bound  estimates  can  be  seen  as  the  sum  of  i.i.d.  random 
terms  which  for  sample  sizes  of  30  or  more  can  be  assumed  normally  distributed  with 
known  (derived  from  the  sampling  process)  variances.  A  95%  confidence  interval  of 
the  obtained  solution  is  computed. 

8.  Computational  Experience 

Computational  results  of  using  Benders  decomposition  and  importance  sampling  for 
two-stage  asset  allocation  problems  can  be  found  in  Infanger  (1991)  [18]  and  Dantzig 
and  Infanger  (1991)  [7]  where  we  report  on  the  solution  of  test  problems  with  up 
to  52  stochastic  parameters  and  a  number  of  universe  scenarios  of  more  than  10^^. 
These  problems  were  formulated  as  two-stage  stochastic  programs.  Using  importance 
sampling  and  sample  sizes  between  200  and  600  very  accurate  results  were  obtained, 
e.g.  the  estimated  95%  confidence  interval  was  less  than  0.8%  on  each  side  based  on 
the  optimal  objective  function  value.  Additional  tests  on  these  examples  showed  that 
the  ratio  of  variance  reduction  obtained  by  using  importance  sampling  versus  crude 
(naive)  Monte  Carlo  sampling  was  about  10“®. 

Inspired  by  these  results  we  implemented  an  earlier  version  of  the  methodology 
described  above  for  the  multi-stage  case  which  did  not  consider  dependency  between 
stages.  Instead  of  the  path  sampling  procedure  for  obtaining  upper  bounds  we  imple¬ 
mented  a  procedure  where  we  sampled  points  rather  than  paths  which  requested  the 
handling  of  an  exponentially  expanding  decision  tree.  Therefore  even  when  we  used 
very  small  sample  sizes,  the  number  of  stages  that  was  practical  to  solve  was  limited. 

We  did  test  up  to  3-stage  problems.  FI3  is  a  3-stage  test  problem  derived  from  a 
2-stage  financial  portfolio  problem  found  in  Mulvey  and  Vladimirou  (1989)  [22].  The 
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problem  is  to  select  a  portfolio  which  maximizes  expected  returns  in  future  periods 
taking  into  account  the  possibility  of  revising  the  portfolio  in  each  period.  There  are 
transaction  costs  and  bounds  on  the  holdings  and  turnovers.  Our  test-problem  covers 
a  planning  horizon  of  3  periods  whereas  the  original  Mulvey-Vladimirou  test-problem 
was  a  2-stage  problem  which  compressed  all  future  periods  into  a  single  second  stage. 
They  solved  the  stochastic  problem  by  restricting  the  number  of  scenarios  in  fl. 

We  assumed  the  returns  of  the  stocks  in  the  future  periods  to  be  independent 
stochastic  parameters  with  3  outcomes  each.  With  13  assets  with  uncertain  returns, 
the  problem  had  26  stochastic  parameters  instead  of  39  because  after  the  last  stage 
decision  was  made,  the  expected  money-value  of  the  portfolio  can  be  evaluated.  The 
number  of  universe  scenarios  was  2.5  - 10*^.  (The  deterministic  equivalent  formulation 
of  the  problem  has  more  than  10'“*  rows  and  a  similar  number  of  columns.)  We 
obtained  an  estimated  optimal  solution  of  the  3-stage  stochastic  problem  using  a 
sample  size  of  only  50  per  stage.  The  optimal  objective  function  value  was  estimated 
to  be  1.10895  with  an  estimated  95%  confidence  interval  of  0.004%  on  the  left  side  and 
0.001%  on  the  right  side  of  the  obtained  objective  function  value.  Thus  the  optimal 
objective  value  lies  within  1.10881  <  2*  <  1.10895  with  95%  probability.  Note  how 
small  the  confidence  interval  is. 

9.  Conclusion 

We  have  demonstrated  how  real-world  multi-period  asset  allocation  problems  can 
be  efficiently  solved  as  multi-stage  stochastic  linear  programs  using  our  approach  of 
combining  Benders  decomposition  and  importance  sampling.  The  numerical  results 
obtained  so  far  are  very  promising:  We  obtained  very  accurate  solutions  for  a  3-stage 
asset  allocation  test-problem  using  remarkably  small  sample  sizes. 
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